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ABSTRACT 

We use cosmological ACDM numerical simulations to model the evolution of the 
substructure population in sixteen dark matter haloes with resolutions of up to seven 
million particles within the virial radius. The combined substructure circular velocity 
distribution function (VDF) for hosts of 10^^ to 10"'^^ M© at redshifts from zero to 
two or higher has a self-similar shape, is independent of host halo mass and redshift, 
and follows the relation dn/dv = (l/8)(vcmax/vcmax,host)~^- Halo to halo variance in 
the VDF is a factor of roughly two to four. At high redshifts, we find preliminary 
evidence for fewer large substructure haloes (subhaloes). Specific angular momenta 
are significantly lower for subhaloes nearer the host halo centre where tidal stripping 
is more effective. The radial distribution of subhaloes is marginally consistent with the 
mass profile for r ^ 0.3rvir, where the possibility of artificial numerical disruption of 
subhaloes can be most reliably excluded by our convergence study, although a subhalo 
distribution that is shallower than the mass profile is favoured. Subhalo masses but 
not circular velocities decrease toward the host centre. Subhalo velocity dispersions 
hint at a positive velocity bias at small radii. There is a weak bias toward more circular 
orbits at lower redshift, especially at small radii. We additionally model a cluster in 
several power law cosmologies of P oc fc", and demonstrate that a steeper spectral 
index, n, results in significantly less substructure. 

Key v^rords: galaxies: haloes - galaxies: formation - methods: N-body simulations - 
cosmology: theory - cosmology: dark matter 



1 INTRODUCTION 

A critical test of the ACDM model is its ability to ac- 
curately predict the evolution of the distribution of "sub- 
haloes" within dark matter haloes, or haloes within haloes. 
The hierarchical formation process of CDM haloes by mul- 
tiple mergers (White & Rees 1978) leaves behind tidally- 
stripped merger remnants that survive as bound subhaloes 
within larger haloes (Ghigna et al. 1998). Subhaloes serve 
as hosts for visible galaxies within clusters, groups, or larger 
galaxies, and so provide a powerful and observable cosmolog- 
ical probe. In cases where dark matter subhaloes may have 
no luminous counterparts, the substructure population can 
be inferred from gravitational lensing studies (e.g. Mao & 
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Schneider 1998; Metcalf & Madau 2001; Chiba 2002; Dalai 
& Kochanek 2002; Mao et al. 2004). 

CDM models predict the abundance of substructure to 
be roughly self-similar, independent of halo mass (Klypin et 
al. 1999; Moore et al. 1999), and to follow a poisson distri- 
bution (Kravtsov et al. 2004a). Subhalo numbers predicted 
by the ACDM model are reasonably matched by observa- 
tions of clusters (Ghigna et al. 2000; Springel et al. 2001; 
Desai et al. 2004; see however Diemand, Moore & Stadel 
2004b). However, observations measure roughly an order of 
magnitude fewer subhaloes in Galactic haloes than in clus- 
ters (Klypin et al. 1999; Moore et al. 1999). Thus, mass is 
either more smoothly distributed on small scales than pre- 
dicted by ACDM cosmology, or Galactic dark matter sub- 
haloes are poorly traced by stars. An element of uncertainty 
in the comparisons with ACDM model predictions is the 
possibility of significant halo to halo variation in the subhalo 
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population that could depend on host mass, merging history, 
or environment. Results from the handfuU of high-resolution 
simulations on galactic scales to date suggest that such vari- 
ance in substructure numbers are significant but much too 
small to account for the apparent discrepancy in Galactic 
subhaloes (e.g. Moore et al. 1999; Klypin et al. 1999; Font 
et al. 2001). Even if observed Galactic dwarfs reside in sub- 
haloes of significantly deeper potential than inferred from 
their stellar velocity dispersions and radial extent, which 
could allow the most massive satellites to match predictions 
(Stoehr et al. 2002; Hayashi et al. 2003; however see also 
Kazantzidis et al. 2004 and Willman et al. 2004), then the 
non-detection of the vast majority of small subhaloes re- 
mains an unsolved problem. However, semi-analytic work 
suggests that baryonic physics causes small haloes to remain 
starless, indicating that observations may be consistent with 
the ACDM model (e.g. Bullock, Kravtsov & Weinberg 2000; 
Benson et al. 2002ab; Somerville 2002). 

By analysing substructure in a large number of dark 
matter haloes, we can measure the range of halo to halo 
variation, and better constrain the uncertainty in the ACDM 
subhalo distribution with the improved statistics. With bet- 
ter numerical resolution and a broad range in host halo 
masses, substructure can be used to place cosmological con- 
straints at new mass scales. Furthermore, detailed dark mat- 
ter simulations provide a theoretical ACDM baseline to link 
subhalo properties with observable galaxy characteristics. 
The conditional luminosity function, which describes the 
number of galaxies of luminosity L ± dL/2 in hosts of a 
given mass, can be combined with simulated subhalo popu- 
lations in order to associate subhalo properties with observ- 
able characteristics (Yang, Mo & van den Bosch 2003; van 
den Bosch, Yang & Mo 2003; see further Vale & Ostriker 
2004) . This allows one to investigate the impact of baryonic 
physics on the galaxy distribution. 

Our simulation set is more sensitive to possible depen- 
dence of the subhalo population on host halo mass or red- 
shift than previous works. Consequently, we can improve 
constraints on whether substructure properties are a func- 
tion of M/Mt, where M* is the characteristic mass of col- 
lapsing haloes defined by the scale at which the rms linear 
density fiuctuation equals the threshold for non-linear col- 
lapse. One might expect that because low-mass haloes were 
mostly assembled much earlier when the universe was more 
dense, are at smaller M/M*, and lie where the power spec- 
trum of mass fluctuations is steeper, that they might have 
significantly different subhalo distributions than more mas- 
sive haloes. If subhaloes within low-mass and low-redshift 
haloes have a higher characteristic density (see e.g. Reed 
et al. 2005 and references therein), then they should be less 
subject to tidal stripping and disruption. Additionally, if the 
infall rate on to the host halo (i.e. merger rate) is different 
from the rate of subhalo destruction, then subhalo numbers 
will evolve with redshift. 

The angular momentum of a cosmic (sub)halo is cru- 
cial to determine the radial distribution of its eventual stel- 
lar component (e.g. Mo, Mao & White 1998; Verde, Oh & 
Jimenez 2002, Van den Bosch et al. 2002) and its collapse 
factor compared to the parent dark matter halo (Stoehr et 
al. 2002). If high angular momentum material is systemati- 
cally stripped from subhaloes this would further inhibit the 
formation of stellar discs in dense environments consistent 



with observations (e.g. Goto et al. 2003). Also, a high col- 
lapse factor of the baryonic component would make the ve- 
locity dispersions of the stars eventually formed lower than 
that associated with the parent halo. This would exacerbate 
the apparent difference between the observed abundance of 
galaxy satellites and the predicted abundance of their host 
haloes, as discussed in detail in Stoehr et al. (2002) and 
Kazantzidis et al. (2004). 

The angular momentum of the dark matter component 
of a subhalo will be relevant to its stellar properties only 
if star formation continues after infall and tidal stripping, 
which requires that sufficient baryons remain in the sub- 
halo. While one generally expects much of a subhalo's gas 
to be stripped upon infall, some of that gas may also be able 
to cool and form stars. Complex processes such as starbursts 
after infall and redistribution of baryonic angular momen- 
tum under tidal infiuences may affect the final spin of the 
baryonic component. Simulations predict starbursts will oc- 
cur in gas-rich satellites as they are tidally stirred by the 
central potential (e.g. Mayer et al. 2001), which may ex- 
plain the continuing star formation in Galactic satellites. 
However, a better understanding of post-infall star forma- 
tion is needed to determine whether the angular momentum 
of the stellar component is correlated with that of the dark 
matter component. 

Subhaloes are particularly sensitive to resolution issues 
(e.g. Moore, Katz & Lake 1996). Dark matter haloes, and 
by extension, subhaloes, have densities that continually in- 
crease toward the halo centre, and so should be very diffi- 
cult to disrupt completely unless numerical discreteness ef- 
fects artificially lower the central density. A subhalo with 
a numerically softened cusp is more easily disrupted by the 
global tidal forces and interactions with other subhaloes that 
strip away the outer regions. Poor spatial resolution and two 
body relaxation lowers central densities, and so may lead to 
subhalo destruction, especially near host halo centres. Sim- 
ulated clusters may generally suffer more from resolution 
issues than galaxies because of their later formation epoch, 
which means that cluster particles will have spent more time 
in low mass haloes (Diemand et al. 2004a). Also, subhaloes 
with highly eccentric orbits are more likely to be disrupted 
since they pass near to the central potential. 

In this work, we present the results of substructure anal- 
yses of 16 ACDM simulated haloes covering three decades in 
mass, from dwarfs to clusters, each with roughly a million 
particles. Our sample includes ten clusters extracted from 
one cosmological volume (CUBEHI) to study cosmological 
variance, and also includes a seven million particle group 
and a four million particle cluster. Some of our haloes are 
well-resolved to redshifts of three or higher, allowing inves- 
tigation of mass or redshift-dependent trends. Additionally, 
we have modelled a cluster in power law cosmologies where 
Pock" to analyse the dependence of the subhalo distribution 
on spectral index, n. 



2 NUMERICAL TECHNIQUES 

2.1 The Simulations 

We use the parallel K-D (balanced binary) Tree (Bentley 
1975) gravity solver PKDGRAV (Stadel 2001; see also Wad- 
sley, Stadel & Quinn 2004) to model sixteen dark matter 
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haloes, further described in Reed et al. (2005); see Table 
1. The CUBEHI run consists of a cube of 432^ particles 
of uniform resolution. The six "renormalized volume" runs 
(e.g. Katz & White 1993; Ghigna et al. 1998) consist of a 
single halo in a high-resolution region nested within a lower 
resolution cosmological volume. Initial conditions for these 
high resolution halo runs are created by first simulating a low 
resolution cosmological volume. Next, a halo of interest is 
identified. To minimise sampling bias, volume-renormalized 
haloes are selected by mass with the only additional con- 
straint that they not lie within close proximity (2 — 3rvir) to 
a halo of similar or larger mass. Then, the initial conditions 
routine is run again to add small-scale power to a region 
made up of high resolution particles that end up within ap- 
proximately two virial radii of the halo centre, while pre- 
serving the original large scale random waves. This process 
is iterated in mass resolution increments of a factor of eight 
until the desired resolution is achieved. We have verified that 
the high-resolution haloes are free from significant contam- 
ination by massive particles. 

Our largest halo has seven million particles and most 
have ~ 10^ particles within the virial radius at redshift zero. 
Due to the high sensitivity of the subhalo distribution to 
numerical resolution effects, we only consider haloes with 
at least 3.5x10^ particles. We adopt a ACDM cosmology 
with Qjn = 0.3 and A = 0.7. The initial density power spec- 
trum is normalized such that erg extrapolated to redshift 
zero is 1.0, consistent with both the cluster abundance (see 
e.g. Eke, Cole & Frenk 1996 and references therein) and the 
WMAP normalization (e.g. Bennett et al. 2003; Spergel et 
al. 2003). We use a Hubble constant of h =0.7, in units of 
100 km s~^ Mpc~^, and assume no tilt (i.e. a primordial 
spectral index of 1). To set the initial conditions, we use the 
Bardeen et al. (1986) transfer function with F = Qm x h. For 
the volume-renormalized runs, we list the effective particle 
number of the highest resolution region rather than the ac- 
tual particle number in Table 1. Numerical parameters are 
consistent with empirical studies (e.g. Moore et al. 1998; 
Stadel et al. 2001; Power et al. 2003; Reed et al. 2003). 

2.2 The Analysis 

The host haloes are defined using the SO algorithm 
(Lacey & Cole 1994) with a tophat overdensity based 
on Eke et al. (1996) where the ACDM virial overdensity, 
Avir, in units of critical density is approximately 100. 
Bound subhaloes are identified using SKID (Stadel 2001; 
[http://www-hpcc.astro.washington.edu/tools/skid.html I , 
which uses local density maxima to identify bound mass 
concentrations independently of environment. SKID iter- 
atively "slides" particles toward higher densities until a 
converged group of particles is found. The radial extent of 
each SKID halo is determined by the distribution of bound 
particles, and no predetermined subhalo shape is imposed. 
The maximum circular velocity of each subhalo, Hcmax, is 
calculated from the measured peak of the rotation curve 
Vc(r) = (GM(< r)/r)°'^. For virialized host haloes, the ratio 
of Vcmax,host/vc,host(rvir) depcuds primarily on M/M^,, and 
is readily estimated from the halo density profile given in 
Reed et al. (2005). Typical values of Vcmax,host/vc,host(rvir) 
are 1.5-1.6 for our dwarfs and 1.05-1.25 for our clusters at 
redshift zero, declining to smaller values at higher redshift 



for each halo. To test for self-similarity between haloes, we 
thus normalize the velocity distribution function (VDF) to 
the peak circular velocities of hosts. We note that a subhalo 
bound to another subhalo will sometimes be catalogued 
as a separate SKID subhalo, particularly when the chosen 
linking parameter r is small. Our tests indicate that the 
velocity distribution function is insensitive to r except 
for the few largest subhaloes. We set r = 2rsoft for the 
CUBEHI run, and r = 4rsoft for all other simulations. We 
consider only subhaloes of 32 or more particles in our 
analyses, consistent with Ghigna et al. (1998) and Springel 
et al. (2001). For analyses of subhalo angular momenta, we 
adopt an empirically motivated minimum of 144 particles 
per subhalo; see section 3.2. Difi'erential plots of binned 
quantities use the median value of each bin, and bin sizes 
are variable in increments of Alogio =0.1, but increased 
when necessary so that no bins are empty. 



2.3 A Resolution Study and the Radial 
Distribution of Subhaloes 

For our volume renormalized runs, GALl, GRPl, and CLl, 
we have analysed 3, 2, and 1 identical lower resolution ver- 
sions, respectively, where the mass resolution is varied by 
increments of 8 x . To find the minimum subhalo circular ve- 
locity down to which our results are complete, Vc.iim, we plot 
the subhalo circular velocity distribution function (VDF) 
for each simulation, and identify the Vcmax below which the 
VDF slope begins to flatten due to incompleteness, as de- 
scribed in Ghigna et al. (2000). In Fig.Q we plot the resolu- 
tion criteria for halo GRPl, marking our conservative Vc,iim 
completeness limits. The agreement of the lower resolution 
versions shows that this technique is sound. Less conserva- 
tive - but still apparently sound - completeness limits for 

1 /3 

our halo sample would have yielded roughly Vc,iim oc N j^^j^ 
for each of our three haloes with multiple resolutions. Nu- 
merical effects on the subhalo population are likely mani- 
fested more strongly for low-redshift subhaloes, since they 
have been subject to potentially disrupting events for more 
time. This implies that measuring Vcmax at z = should still 
be valid at higher redshift. We have verified that the z = 
Vc,iim for halo GRPl remains valid at z = 1 where the high- 
est resolution version still has 3.5 x 10^ particles. Addition- 
ally, we have considered whether the completeness limits, 
which were obtained using subhaloes within the entire virial 
volume, may change at small radii, where numerical prob- 
lems are expected to be stronger. We find that even when 
the sample is limited to r < 0.3rvir, the Vc,iim derived from 
the entire virial volume as described above, remains valid. 

The radial distribution of the subhalo population, 
shown for GRPl in Fig.|21 has a slope equal to or shallower 
than the slope of the density profile for each resolution at 
all radii. At very small radii, subhaloes are highly deficient, 
possibly due at least in part to artificial disruption by the 
tidal forces from the central mass concentration. When a 
subhalo migrates inward via dynamical friction to a radius 
where its central density is lower than the local density of 
the host halo, it will likely be disrupted (e.g. Syer & White 
1998). Numerical disruption enhanced by the strong tides 
near host centres could affect subhaloes. The flattened slope 
of the radial subhalo distribution of GRPl subhaloes rela- 
tive to the mass profile interior to roughly 0.2rvir for our 
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Table 1. Summary of our halo sample at redshift zero. For volume renormalized runs, the mass (h^'^MQ) and particle number of 
the central halo is listed. Vcmax,host is the host peak circular velocity (km s~^). Np eff is the effective particle number based on the 

^clim is i'h'2 subhalo peak circular velocity (km s~^) above which numerical disruption is 
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Figure 1. VDF for halo GRPl (at z = 0) with 7 million parti- 
cles, and for 2 lower resolutions, with particle mass incrementally 
increasing by a factor of eight. The circular velocity completeness 
limit Vcjini is marked for each resolution with a vertical line. 



highest resolution, and interior to O.Srvir for our lowest reso- 
lution implies that disruption and/or stripping of subhaloes 
is important in the halo central region. The increase in ra- 
dius of the break in subhalo slope with decreasing resolution 
suggests that numerical disruption, if present, is worse for 
lower resolutions. However, we caution that the location of 
the break is not well defined because of poisson uncertain- 
ties. Given the uncertainties, it is not possible to reliably 
separate spurious from real disruption that may be present 
in our simulations at small radii. Thus, we have no evidence 
that the central substructure number density has converged 
with resolution. Increasing the mass resolution by a factor of 
eight results in roughly a factor of ~2-2.5 more subhaloes 
at a given radius beyond roughly O.Srvir, though there is 
substantial noise in this estimate. At larger radii, where nu- 
merical effects are less important, a subhalo radial distribu- 



Figure 2. The z = radial distribution of subhaloes in GRPl 
for each of the 3 different resolutions (solid, long-dashed, short- 
dashed, are the highest, mid, and lowest resolution runs, respec- 
tively). Thin lines have the Vcjim resolution limit of the lowest 
resolution run applied, la poisson error bars are shown. The dot- 
ted line is the particle distribution with arbitrary normalization. 



tion that is shallower than the mass profile is favoured in 
our data, and is also reported by Diemand et al. (2004b), 
Gill, Knebe & Gibson (2004a), Stoehr et al. (2003), Gao et 
al. (2004ab), De Lucia et al. (2003), and Nagai & Kravtsov 
(2005). However, a radial subhalo slope equal to the mass 
profile slope is not ruled out except in the central region 
where numerical effects may dominate. Note that if we use 
only a minimum mass cut without imposing the circular 
velocity completeness limit, a more significant anti-bias in 
the radial distribution is found (see also Nagai & Kravtsov 
(2005) 
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3 RESULTS 

3.1 The Circular Velocity Distribution Function 

(VDF) 

In Fig. 1^1 we plot the cumulative VDF for all haloes at red- 
shift zero. Our completeness and resolution limits exclude 
the majority of subhaloes identified by SKID, leaving of the 
order of 100 subhaloes in most hosts. Fig. 2] shows the dif- 
ferential VDF for our entire set of subhaloes at ^z^2 for 
hosts of more than 3.5 x 10^ particles. The subhalo VDF is 
well-approximated by a power law with slope and normal- 
ization given by 



dn/dv=-( 



o Vcniax.host 



(1) 



over the range of approximately 0.07vcmax,host to 
0.4vcmax,host with halo to halo scatter of a factor of 
roughly two to four. When we individually consider the 
VDF for hosts of different mass, we find no evidence of 
mass dependence on the VDF amplitude or slope in our 
data, and similarly we detect evidence of weak or zero 
redshift dependence (see below). This implies that the 
number density of subhaloes is approximately self-similar, 
independent of the mass and redshift of the host halo. Here 
we caution that a larger halo sample would be needed to 
identify any weak trends (e.g. an increase in substructure 
abundance with halo mass, reported by Gao et al. 2004a) 
as they would be masked by the large halo to halo scatter. 
The farthest low VDF outlier halo is the largest cluster in 
the CUBEHI simulation a.t z = 1, which has approximately 
4.5 X 10^ particles at the time. In Fig. |S] we present the 
VDF of our highest resolution halo to redshift z ~ 4. There 
is little or no evolution in VDF slope or normalization, 
although there are fewer large subhaloes at high redshift, 
as seen by the lack of data points in the VDF beyond large 
Vcmax/vcmax,host iu high-redshift hosts. Thus, large subaloes 
are deficient until lower redshifts when they either infall on 
to the host or are formed as merger products of existing 
subhaloes. 

In Fig. 1^1 we present the non-normalized differential 
VDF for our entire set of subhaloes, which displays little or 
no redshift evolution in comoving space. The subhalo VDF 
is fit by 



dn/dv = 1.5 X 10V;rmax(li^Mpc"^km~^s) 



(2) 



for subhaloes of 20-300 kms~^ in hosts of lO" to 10^" 
Mq and redshifts from zero to two, with scatter of a fac- 
tor of ~ 2 — 4. This slope of this non-normalized VDF 
is steeper because lower mass galaxies in the sample have 
a larger ratio of Vcmax,host/vc,host(rvir). The dependence 
on Vcmax. host /vc, host (ivir) implies that this non-normalized 
VDF is unlikely to hold universally. Fig. |7| shows that the 
non-normalized subhalo VDF for GRPl has weak or no evo- 
lution to 2: ~ 4. 

We also calculate the VDF of friends-of-friends (FOF, 
Press & Davis 1982; Davis et al. 1985) haloes, from which 
the population of subhaloes must have descended. The FOF 
algorithm identifies virialized haloes that are not part of 
larger virialized objects. In Fig. |H| we plot the VDF for all 
FOF haloes in the CUBEHI simulation adopting a linking 
length of 0.2 times the mean particle separation. We have 
normalized the FOF VDF by multiplying it by the mean 
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Figure 3. The 2 = cumulative subhalo VDF with the complete- 
ness limits V(, ii,ji applied, and considering only subhaloes with 
32 or more particles, plotted as a function of Vcmax/vcmax host- 
Solid lines are renormalized volumes (green, light blue, magenta, 
red, and blue for CLl, GRPl, GRP2, GALl, and DWFl, respec- 
tively). Dashed lines (black) are the ten clusters in the CUBEHI 
simulation. 
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Figure 4. Normalized subhalo VDF for all haloes with 
Np > 3.5 X 10^ at redshifts 0, 0.5, 1, and 2, normalized to 
Vcmax,host s-ud a virial volume of unity. Solid lines (colours as in 
Fig.|2J are renormalized volumes, dashed lines (black) are the ten 
clusters in the CUBEHI simulation. Heavy solid line corresponds 

to dn/dv = ^(vcmax/vcmax.host)"'*- 
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Figure 5. Evolution of the normalized subhalo VDF for halo 
GRPl. Subhalo cireular velocities are normalized to Vcmax host 
and a virial volume of unity. 



virial overdensity, Pviriai/p, at each redshift. Differences be- 
tween the subhalo and halo VDFs are generally smaller than 
a factor of three until z ^ 7, in agreement with z = results 
by Gao et al. (2004a). At z > 5, the FOF VDF drops rapidly 
with increasing redshift, presumably due to the fact that it 
samples the steep drop-off regime of the mass function at 
high redshift. Higher redshift simulations are needed to test 



whether the lower FOF VDF at 



7 results in small sub- 



halo numbers for extremely high-redshift hosts. 



3.2 Subhalo Angular Momenta 

We analyse the subhalo angular momenta given by the spin 
parameter, A = LE^''^ /GM^^'^ (Peebles 1969), for sub- 
haloes at redshift zero in our two highest resolution simula- 
tions, GRPl and CLl. Angular momentum L is calculated 
with respect to the centre of mass of the subhalo. Energy 
E = Bkinetic + -Bpotentiai is summcd ovcr all subhalo par- 
ticles. Each central subhalo and subhaloes with more than 
50,000 particles (for computational speed) are excluded. In 
our simulations, median subhalo A increases with decreas- 
ing particle number once below ~100 particles, indicating an 
upward bias for poorly resolved haloes. We thus limit our 
spin analyses to subhaloes containing 144 or more particles, 
which results in no dependence of median spin on subhalo 
mass. Several studies have found the lognormal function to 
be a good description of halo spins 



p(A)dA : 
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ctaVStt 
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ln2(A/Ao)\ d-^ 
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(3) 



(Barnes & Efstathiou 1987; Ryden 1988; Cole & Lacey 1996; 
Warren et al. 1992; Gardner 2001; van den Bosch et al. 2002; 
Vitvitska et al. 2002; Cohn et al. 2004; Peirani, Mohayaee 
& Pacheco 2004; Aubert, Pichon & Colombi 2004)) where 



Figure 6. Differential subhalo VDF for all haloes with Np > 
3.5 X 10^ at redshifts 0, 0.5, 1, and 2. Solid lines (colours as in 
Fig. 121 are renormalized volumes, dashed lines (black) are the ten 
clusters in the CUBEHI simulation. The heavy solid line denotes 
dn/dv = 0.5 X 10*V;r„jax (Eq.|2j. Volume units are comoving. 



Ao and g\ axe fit parameters. Fig. [^compares the histogram 
of subhalo spins with the best-fitting x^ lognormal function 
given by (Ao,aA) = (0.0235, 0.54) and (0.0238, 0.73) for 
halo GRPl and CLl, respectively. Median and arithmetic 



; = 0.027 for GRPl 
0.026 for CLl sub- 



average values of Amcd = 0.024 and Aa 
subhaloes, and Amcd = 0.024 and Aavg 
haloes. Subhalo spins are significantly smaller than spins 
of a sample containing 1.5 x 10* field haloes selected from 
the CUBEHI volume using the spherical overdensity (SO) 
algorithm, which has Amod = 0.037 and (Ao,o-a) = (0.037, 
0.57) over the mass range of 1.9 x 10^° - 6.5 x 10^^ (144- 
50000) particles. Within this field halo sample and within 
the subhalo samples, spins have no significant mass depen- 
dence, so comparisons between differing mass scales of dif- 
ferent simulations should be valid. Our field angular mo- 
menta are consistent with a number of recent studies that 
find Ao = 0.035-0.046 for virialized ACDM haloes (Bullock 
et al. 2001; van den Bosch et al. 2002; Vitvitska et al. 2002; 
Colin et al. 2004; Peirani et al. 2004; Aubert et al. 2004). 
Note that some of these studies use a slightly different def- 
inition of spin introduced by Bullock et al. (2001) that has 
some dependence on the density profile. 

Given that our subhaloes are in high-density environ- 
ments, it is likely that their spins are lowered by the removal 
of high angular momentum material, which should be most 
vulnerable to tidal stripping. We verify that stripping of 
outer material has the potential to lower spins by the re- 
quired amounts by measuring spins for the central regions 
of SO haloes in the CUBEHI volume. Here we find that the 
central SO region containing 15 x the usual virial overden- 
sity has Amed = 0.025, which is similar to Amod for subhaloes. 
This overdensity generally contains the central ~ 15-30% of 
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Figure 7. Subhalo VDF in comoving coordinates for the group 
GRPl, which still has over half a million particles at z = 3.8. 1 
a poisson error bars are shown. 



Figure 8. VDF for FOF haloes in the CUBEHI simulation plot- 
ted for ^ 2 ^ 15. Heavy solid line is the fit to the subhalo pop- 
ulation given by Eq. 2 (non- normalized). For reference, we nor- 
malize the FOF VDF to the virial overdensity as a function of 
redshift by multiplying by a factor of Pvirial/p- 



the SO mass. We have also determined that the SO haloes 
whose central material have the lowest spins relative to the 
surrounding halo material are the high-spin SO haloes that 
comprise the tail of the lognormal distribution. This sug- 
gests that subhaloes with initially high spins are likely to 
be the largest contributors to the removal of subhalo angu- 
lar momentum. Interestingly, the central regions of low-spin 
SO haloes (A << Amed) generally have higher spin than the 
surrounding halo material. 

An alternative explanation for the low spins of sub- 
haloes may be related to the fact that (sub) halo mass growth 
is generally halted upon infall into a larger halo. For the 
case of field haloes, spins decrease over time in the absence 
of major mergers, particularly in the period immediately 
following a major merger event (e.g. Vitvitska et al. 2002; 
Peirani et al. 2004). If a similar decrease in spin occurs in 
the high density environments within haloes, then it could 
contribute to the low subhalo spins. More study is needed 
to determine whether this mechanism for lowering spin is 
effective for subhaloes. 

The radial dependence of subhalo spins shown in Fig. 1101 
shows that local environment directly affects subhalo spin. 
Median subhalo spin decreases with decreasing radius from 
A ~ 0.025 - 0.03 near rvir to A ~ 0.02 at O.Srvir, where 
tidal effects are greater. The transition in spins from sub- 
haloes to field haloes at rvir is smooth, and spins of SKID 
haloes found between rvir and 2rvir have weak or no fur- 
ther radial dependence. Note that the median SKID halo 
spin Amed — 0.03 between rvir and 2rvir is somewhat smaller 
than for SO haloes selected from the uniform volume. We 
have tested that when the field SO haloes of the large vol- 
ume (CUBEHI) are instead selected with SKID, their spins 
are smaller than with SO, but only by 5 — 10%, which is not 
enough to explain the small spins found in the high-density 



regions. This suggests that field haloes in the high-density 
environments just outside of larger virialized haloes have 
spins that are ~10% lower than for the global population. 
Spins in high-density regions could be reduced due to a con- 
tribution from subhaloes whose orbits have previously taken 
them within the virial radius (Gill, Knebe & Gibson 2005), 
causing stripping of high angular momentum material. Ad- 
ditionally, tidal interactions with the massive neighbouring 
halo (Gnedin 2003; Kravtsov, Gnedin & Klypin 2004b) may 
also have some impact on spin. This relatively weak envi- 
ronmental dependence on field halo spins is still consistent 
with Lemson & Kauffmann (1999), who found no difference 
in spins of virialized haloes in mild over(under)-densities. 

The radial trend in spins could be observable if reliable 
spins can be estimated for a large sample of cluster galax- 
ies or satellites in external systems, and if spin of baryonic 
material can be used as a tracer of subhalo spin. Galaxies 
that form from low spin material should have larger collapse 
factors. If star formation occurs after subhalo angular mo- 
mentum has been lowered, then galaxies nearer to host cen- 
tres may thus have smaller radial extent than galaxies near 
or outside the virial radius, although complex baryonic pro- 
cesses related to star formation may dominate over any po- 
tential spin-induced trend in stellar distribution (see Stoehr 
et al. 2002; Hayashi et al. 2003; Kazantzidis et al. 2004). 
The central galaxy would likely deviate from this trend since 
it undergoes late mergers, which have been shown to in- 
crease halo angular momentum (Gardner 2001; Vitvitska et 
al. 2002). 

A radial trend in subhalo spins is unlikely to be eas- 
ily observable unless the stellar component of a galaxy is 
assembled subsequent to infall and stripping of high spin 
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Figure 9. Histograms siiow the z = distribution of spin subhalo 
parameters within halo GRPl (solid) and halo CLl (dashed) for 
subhaloes with 144 or more particles. Spins of SO field haloes 
from the CUBEHI volume are shown by dotted lines. Curves are 
fits to the lognormal distribution described in text. 



material. If gas infall, star formation, and disc growth all 
cease soon after infall, then spins of galaxies hosted by sub- 
haloes would only reflect the spin of the original dark matter 
halo. If however, sufficient gas remains after accretion into 
another dark matter halo that there is some continuing star 
formation (as for Galactic satellites), then the stellar and 
gaseous distribution may have a smaller spin than that of 
the dark matter host prior to infall. Numerical simulations 
are needed to determine the extent that subhalo baryonic 
spins reflect the spins of their host dark matter subhaloes. 



3.3 Subhalo Mass Function 

We plot the subhalo mass function at z = (Fig. IIH and 
aX z = 1 (Fig. I12II . The steep drop on the low mass end for 
each halo is due to exclusion of haloes with Vcmax < Vc,iim 
for each simulation. The subhalo mass function is indepen- 
dent of halo mass, a result also seen in De Lucia et al. 2004, 
although a weak mass-dependence such as that reported by 
Gao et al. 2004a is not ruled out by our data. The Sheth and 
Tormen (1999) function normalized by a factor equal to the 
virial overdensity, is plotted for reference. The Sheth and 
Tormen function, a modification of the Press and Schecter 
(1974) formalism, is an excellent match to the CUBEHI FOF 
mass function at low redshifts (Reed et al. 2003) . The factor 
of approximately two to three offset between the subhalo 
mass function and the FOF mass function is independent 
of meiss and redshift, which implies that the stripping effi- 
ciency of subhaloes is largely mass and redshift independent, 
though infall timing and evolution of the global mass func- 
tion may also affect the subhalo mass function. 




Figure 10. Radial dependence of spin parameter at z = 0. Solid 
(dashed) histogram is median of subhalo spins for host GRPl 
(CLl) 



3.4 Radial Distribution and Subhalo Orbits 

In Fig. 1181 we show the distribution of subhalo Vcmax versus 
radial position for a range of redshifts. There are no clear 
radial trends in Vcmax for z ^2. However, Fig. ll4l shows that 
subhaloes near the centres of their hosts tend to have lower 
median masses than subhaloes at larger radii, a trend con- 
sistent with e.g. Taffoni et al. (2003) (see also De Lucia et 
al. 2004; Taylor & Babul 2004; Kravtsov et al. 2004b; Gao 
et al. 2004a), indicative of tidal stripping near halo centres 
(e.g. Tormen, Diaferio & Syer 1998). This suggests that the 
morphology-radius relation seen in clusters (e.g. Whitmore 
& Gilmore 1991) cannot be explained by a Vcmax-radius re- 
lation. We note that the region where the median mass is 
smallest (less than ~0.3rvir is also the least numerically ro- 
bust, discussed in section 2.3. 

In Fig. 1151 we show the circularity of subhalo orbits, 
L/Lcirc, which is the angular momentum that a subhalo 
with a given orbital energy would have if it were on a cir- 
cular orbit. Orbits are calculated from a snapshot position 
and velocity of each substructure applied to a static spheri- 
cal approximation of the host halo potential, which is com- 
puted from the mass profile as in Ghigna et al. (1998). Our 
subhaloes have mean circularities between 0.6 to 0.7 where 
sampling is high, with a redshift zero mean L/Lcirc = 0.64 
for r < rvir- The subhalo circularity increases weakly at 
small radii for lower redshifts, suggesting circularisation of 
orbits over time as seen in simulations by (Gill et al. 2004). 
Also, there are a larger number of subhaloes in the nearly 
circular orbits than in the most radial orbits, especially at 
low redshift. Both of these trends may be a result of dis- 
ruption or heavy stripping of subhaloes on extremely radial 
orbits since they pass nearer the cluster central potential. 
This effect would be greatest for subhaloes with small apoc- 
entres. In Fig. 1161 we present the circularity as a function of 
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Figure 11. The mass function of the subhalo population for each 
of our haloes at z = 0. Heavy dark line is the Sheth and Tormen 
(1999) prediction for virialized haloes, which closely matches the 
low redshift FOF mass function, and is normalized by a simple 
factor equal to the virial overdensity (pvirial/p) for reference. 




I . I I I I 
10 11 12 

Logi„{M/h->Mj 



Figure 12. Subhalo mass function, as in Fig. 1111 except at z = 1. 



Vcmax/vcmax,host. Again, there are no strong trends, except 
that small subhaloes with very radial orbits appear to be 
relatively deficient, especially at low redshift. 

Fig. 1171 shows the three dimensional velocity disper- 
sion of subhaloes, ctsd, as a function of radius at red- 
shift 0. The velocity dispersion is flat or slowly decreasing 
with radius at ^ 0.5rvir. At smaller radii, ctsd generally in- 
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Figure 13. The Vcmax/vcmax.host distribution of subhaloes plot- 
ted against radial position at redshifts (top left), 0.5 (top right), 
1 (bottom left), and 2 (bottom right). Individual histograms rep- 
resent the median Vcmax/vcmax.host from each simulation with 
colours as in Fig. U] (e.g. dashed (black) histogram is an average 
of 10 CUBEHI clusters). Coloured points are from the normalized 
runs with same colours as before, and black points are from the 
CUBEHI run. The subhalo associated with the potential center 
is excluded from each halo. Subhaloes from hosts with less than 
3.5 X 10^ particles are again excluded. 
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Figure 14. The M/Mi^Qst distribution of subhaloes plotted 
against radial position at redshifts (top left), 0.5 (top right), 
1 (bottom left), and 2 (bottom right). Individual histograms rep- 
resent the median M/Mjjost from each simulation as in Fig. 1131 
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Figure 15. The circularity of subhalo orbits, L/Ldrc vs. ra/rvir, 
for the same haloes and redshifts as in Fig. I13II14I L^rc is the 
orbital angular momentum of a subhalo on a circular orbit of 
a given orbital energy, ra is the subhalo apocentre. Histogram 
represents the average L/L^irc- 
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Figure 16. The circularity of subhalo orbits, L/Ldrc vs. 
(vcmax/vcmax host)) for th^ Same redshifts as in Fig. 1131^151 here 
limited to subhaloes within r^ir- Average circularity given by his- 
togram. 
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Figure 17. Top: 0-30 vs. radius at z = for each halo. Bot- 
tom: Subhalo velocity bias, b = Tsub.SD/o'SD.dm- Line types and 
colours are the same as Fig.O 



creases moderately toward the centre, reaching less than 1.5- 
2 cr3D(rvir). Subhalo orbital motions suggest a velocity bias 
b=(7sub,3D/o"3D,dm (bottom panel of Fig. I17II that increases 
with decreasing radius, reaching ~ 1.2 by ~ 0.2rvir, simi- 
lar to that found in e.g. Colin, Klypin, & Kravtsov (2000), 
Ghigna et al. (2000), and Diemand et al. (2004b). Diemand 
et al. attribute both their central velocity bias and their 
central spatial anti-bias to tidal destruction of slow moving 
subhaloes near the cluster centre. Note that in a dynamically 
relaxed system, a spatial anti-bias automatically results in a 
positive central velocity bias (see e.g. van den Bosch 2004). 
At radii ^ O.Srvir, there is also a hint of a weak velocity 
anti-bias of ~ 10% in many haloes, although subhalo veloc- 
ities are consistent with no large radius anti-bias when the 
range of halo to halo scatter and uncertainty due to small 
numbers of subhaloes is considered. In Fig. 1181 we plot the 
subhalo CT3D for GRPl to redshift 4 against 0-30 for particles. 
The subhalo velocity bias and a^o are each consistent with 
no redshift evolution. 



3.5 Power Law Cosmologies 

In order to examine the effects of the power spectral slope 
index, n, we have simulated a single renormalized volume 
cluster with a range of values for n. Here the initial density 
fluctuation power spectrum is given by P oc A:", normal- 
ized to og, = 1.0 as in the ACDM simulations. We plot the 
subhalo VDF for the cluster with initial conditions given 
by n =0, -1, -2, and -2.7 in Fig. 1191 There is a clear and 
strong trend that steeper power spectra have less substruc- 
ture. This is a direct result of the fact that the shallower 
spectra have more small scale power relative to large scale 
power than the steeper spectra. Additionally, subhaloes form 
earlier in cosmologies with flatter power spectra, and have 
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Figure 18. Top: Evolution of velocity dispersions for the 7 mil- 
lion particle group GRPl. Thick lines are (T3D of the subhalo 
distribution. Thin lines are (T30 for particle distribution. Bottom: 
Subhalo velocity bias for same redshifts. 
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Figure 19. Subhalo VDF for our cluster with initial power spec- 
trum given by P oc fc". Heavy solid line is a fit to ACDM haloes. 



higher characteristic densities and steeper density profiles 
(e.g. Reed et al. 2005), making them less vulnerable to dis- 
ruption. 



4 CONCLUSIONS 

(1) A "universal" subhalo VDF independent of host mass 
and redshift: The subhalo follows the self-similar relation: 
dn/dv = (1/8) (vcmax/vcmax, host )""* with a factor of two to 
four scatter. 

(2) Subhalo spins decrease toward the host centre. The 
radial dependence of subhalo spins is likely explained by 
the increased vulnerability of high-spin material to tidal 
stripping and disruption. Galaxies that form from low angu- 
lar momentum subhaloes in high-density regions after some 
stripping has occurred would have a larger baryonic collapse 
factor. If we assume the size of the stellar distribution is 
proportional to A (as in e.g. Kravtsov et al. 2004b), then 
stellar systems could be ~50% smaller if they form in heav- 
ily stripped central subhaloes. This would require that star 
formation continue after the subhalo has been stripped of 
high angular momentum material, as discussed previously. 
The circular velocity of the final stellar distribution could be 
lowered by a similar amount if the stars collapse to a small 
radius, assuming a central subhalo density slope near r~^ 
as found in simulations by Kazantzidis et al. (2004). Stoehr 
et al. (2002) proposed that the the apparent deficit of lo- 
cal group satellites could be explained if satellites actually 
reside in large subhaloes (further supported by Hayashi et 
al. 2003) , which would allow the VDF of large Local Group 
satellites to match ACDM predictions. The small radii of the 
stellar distribution would give them low stellar velocity dis- 
persion relative to Vcmax of their dark matter subhalo hosts. 
Baryonic spins of A < 0.01 would be required for the gas 
to collapse to the required ~1% of the pre-stripped virial 
radius of the dark subhalo (Kazantzidis et al. 2004) needed 
to match Local Group satellites with ACDM predictions. 
However, only ~ 5% of our subhaloes have spins of A ^ 0.01, 
and 90% have A > 0.014, so such a solution appears unlikely 
if the steep central substructure density profiles found by 
Kazantzidis et al. are correct. 

(3) Evidence for redshift dependence of the VDF: The 
shape and amplitude of the subhalo VDF has little or no 
trend with mass or redshift. However, fewer subhaloes lie 
at large Vcmax/vcmax,host in high-redshift hosts. This implies 
that dark matter haloes are not populated with large sub- 
haloes until lower redshifts. 

(4) Could the VDF be a function of M/M^ ? Though 
our data is consistent with only a weak or no trend of the 
VDF on host mass and redshift, a trend of less substructure 
in large M/M, haloes would not be surprising, although it 
remains somewhat speculative until simulations are able to 
probe a larger mass and redshift range. The subhalo pop- 
ulation approximately follows some fraction of the univer- 
sal halo population. Thus, at large M/M^, where the mass 
function is steep, the subhalo mas function should also be 
steep. The reason we see no such trend among redshift zero 
haloes in our sample may be because M/K'Lj, for our simu- 
lated haloes is never much larger than unity, below which the 
slope of the field mass function has little mass dependence. 

A dependence of substructure internal density profiles 
on M /h'h could also cause a trend for large M /M^ hosts to 
have decreased substructure. Virialized haloes have density 
profiles that become less centrally concentrated with increas- 
ing M/Mf (see e.g. Reed et al. 2005), so subhalo density 
profiles are likely to have a similar dependence. The lower 



12 Reed et al. 



central densities of low M /M^ subhaloes should make them 
more vulnerable to stripping and major disruption during 
and after virialization, even though their hosts also have 
lower central densities. The potential decrease in substruc- 
ture at large M/M, may be manifested in our highest red- 
shift data (Fig. EJ where M/M, ~ 1000. We note however 
that our power law cosmology shows that the substructure 
abundance is also strongly dependent on the spectral slope. 
The fact that low-mass haloes sample steeper portions of 
the density power spectrum might cancel out any substruc- 
ture dependence on Af/M«, though the range of spectral 
slopes present over our halo mass range is quite small. Also, 
because low mass haloes generally form earlier, one might 
expect them to have less substructure because their sub- 
haloes have been subject to tidal stripping and disruption 
for longer periods of time. 

(5) Is there a central subhalo "anti-bias"? The slope of 
the subhalo number density is marginally consistent with 
that of the mass density for r ^ 0.3 rvir, although a subhalo 
distribution that is shallower than the dark matter profile 
is favoured. More high-resolution simulations are needed to 
conclusively demonstrate a subhalo anti-bias at large radii. 
At smaller radii, subhaloes are deficient, but artificial nu- 
merical effects cannot be excluded. It is not surprising that 
more tidal stripping has occurred for subhaloes at small host 
radii, as reflected in the mass-radius trend of Fig. ll4l and seen 
also by e.g. Gao et al. (2004a). One would also expect that 
the removal of sufficiently large masses would lead to a radial 
trend in Vcmax, although such a radial trend may be masked 
by the preferential destruction of low Vcmax haloes at small 
radii in our results. In any case, since subhaloes are likely to 
have steep cuspy density profiles down to less than 1% rvir 
upon initial infall, as suggested by numerous studies of the 
density profiles of virialized haloes (e.g. Reed et al. 2005), 
then their centres will be highly resistant to tidal disruption. 
Thus, any appearance of subhalo "anti-bias" in current gen- 
eration simulations is likely to be simply a manifestation of 
a radial trend in subhalo mass (or Vcmax). The imposition 
of an arbitrary minimum mass (or minimum Vcmax) as re- 
quired by resolution constraints can give the false impression 
of "missing" central subhaloes (see also Gao et al. 2004b for 
a detailed discussion). 

(6) Subhaloes have no strongly preferred orbits. Orbital 
properties show no strong trends with respect to radius, red- 
shift, or Vcmax. Howcvcr, the subhaloes on the most highly 
eccentric orbits become less abundant over time, likely due 
to disruption or heavy stripping by the central potential of 
the host, but this only affects a small fraction of total sub- 
haloes. Subhaloes with small apocentres are most strongly 
affected. This likely has an artificial numerical cause wherein 
tidally affected subhaloes either lose enough material that 
they drop below resolution constraints, or are completely 
disrupted due to the effective density ceiling imposed upon 
simulated subhaloes. Our results also suggest that subhaloes 
have a positive velocity bias at small radii and little or no 
velocity bias at large radii. It is not clear how the addition 
of baryons in the form of gas and stars would affect orbital 
kinematics and distribution of the subhalo population. 

(7) Cosmological variance is too small to account for 
"missing" satellites. The subhalo VDF has a halo to halo 

cosmological variance of a factor of roughly two to four. 
This means that the apparent problem of overpredicted local 



group satellites cannot be solved by invoking cosmological 
variance. 
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